--------------------------------------------------------------------------------
-- This package is a thick Ada binding to the Multiple Precision Floating Point
-- Reliable Library. See <http://www.mpfr.org/>.
--
-- Copyright (C) 2009 Vincent DIEMUNSCH.
-- GNU General Public License.
--
-- This program is free software: you can redistribute it and/or modify it under
-- the terms of the GNU General Public License as published by the Free Software
-- Foundation, either version 3 of the License, or (at your option) any later
-- version.
--
-- This file is distributed in the hope that it will be useful, but WITHOUT ANY
-- WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-- A PARTICULAR PURPOSE. See the GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License along with
-- this program.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------------
---

with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Text_IO;                               use Ada.Text_IO;
with Ada.Strings.Fixed;
with GMP.Integers.Aliased_Internal_Value;

package body MPFR.Floats is

   use type Specific.mpfr_prec_t;
   use type Specific.mp_exp_t;

   ------------------------
   -- Global Variables : --
   ------------------------

   default_rounding_mode  : mpfr_rnd_t  := GMP_RNDN;
   required_min_precision : mpfr_prec_t := mpfr_prec_t'First;

   -------------------------------------
   -- Precision in Digits and in Bits --
   -------------------------------------

   -- See ARM G.2.2 "Model-Oriented Attributes of Floating Point Types" :
   -- S'Model_Mantissa = ceiling(d * log(10) / log(T'Machine_Radix)) + g
   -- where d is the requested decimal precision of T, and g is 0 if
   -- T'Machine_Radix is a positive power of 10 and 1 otherwise.

   -- T'Machine_Radix is 2 (it is 10 in the case of Binary Coded Decimal),
   -- and therfore g = 1.

   Log2_10 : constant Float := 3.321_928_094_887;

   function Prec_in_Bits (Decimal_Digits : Positive) return mpfr_prec_t is
   begin
      return mpfr_prec_t (Float'Ceiling (Float (Decimal_Digits) * Log2_10)) + 1;
   end Prec_in_Bits;

   function Prec_in_Digits (prec : mpfr_prec_t) return Positive is
   begin
      return Positive (Float'Floor (Float (prec - 1) / Log2_10));
   end Prec_in_Digits;

   pragma Inline (Prec_in_Bits, Prec_in_Digits);

   ---------------------------
   -- Set_Default_Precision --
   ---------------------------

   procedure Set_Default_Precision (Decimal_Digits : Positive) is
   begin
      mpfr_set_default_prec (Prec_in_Bits (Decimal_Digits));
   end Set_Default_Precision;

   ---------------------------
   -- Get_Default_Precision --
   ---------------------------

   function Get_Default_Precision return Positive is
   begin
      return Prec_in_Digits (mpfr_get_default_prec);
   end Get_Default_Precision;

   -------------------------------------
   -- Unchecked_Precision_Adjustement --
   -------------------------------------

   procedure Unchecked_Precision_Adjustement
     (x              : in out MPFR_Float;
      Decimal_Digits : Positive)
   is
   begin
      mpfr_prec_round
        (x.Value,
         Prec_in_Bits (Decimal_Digits),
         default_rounding_mode);
   end Unchecked_Precision_Adjustement;

   procedure Unchecked_Precision_Adjustement
     (x              : in out MPFR_Float;
      Error          : MPFR_Float;
      Decimal_Digits : Positive := 3)
   is
      Prec_in_Decimal_Digits : Positive range 3 .. Positive'Last;
   begin
      if Decimal_Digits < 3 then Prec_in_Decimal_Digits := 3;
                            else Prec_in_Decimal_Digits := Decimal_Digits;
      end if;
      if mpfr_cmpabs (Error.Value, x.Value) < 0 then
         mpfr_prec_round
           (x.Value,
            Prec_in_Bits (Prec_in_Decimal_Digits) +
            mpfr_prec_t'Base (mpfr_get_exp (x.Value) -
                              mpfr_get_exp (Error.Value)),
            default_rounding_mode);
      end if;
   end Unchecked_Precision_Adjustement;

   -------------------------------
   -- Set_Default_Rounding_Mode --
   -------------------------------

   function To_mpfr_rnd_t (Rounding : Rounding_Mode_Kind) return mpfr_rnd_t is
   begin
      case Rounding is
         when Round_to_Nearest            => return GMP_RNDN;
         when Round_Toward_Zero           => return GMP_RNDZ;
         when Round_Toward_Plus_Infinity  => return GMP_RNDU;
         when Round_Toward_Minus_Infinity => return GMP_RNDD;
      end case;
   end To_mpfr_rnd_t;

   procedure Set_Default_Rounding_Mode (Rounding : Rounding_Mode_Kind) is
   begin
      default_rounding_mode := To_mpfr_rnd_t (Rounding);
   end Set_Default_Rounding_Mode;

   ----------------------------------------------------
   -- Set_Precision_Check & Inhibate_Precision_Check --
   ----------------------------------------------------

   procedure Set_Precision_Check (Decimal_Digits : Positive) is
   begin
      required_min_precision := Prec_in_Bits (Decimal_Digits);
   end Set_Precision_Check;

   procedure Inhibate_Precision_Check is
   begin
      required_min_precision := mpfr_prec_t'First;
   end Inhibate_Precision_Check;

   -------------------
   -- To_MPFR_Float --
   -------------------

   function To_MPFR_Float
     (Image          : String;
      Decimal_Digits : Natural := Default_Precision)
      return           MPFR_Float
   is
      -- This implementation is derived from the one of GMP.To_Unbounded_Integers
      use Ada.Strings.Fixed;
      mpf           : mpfr_t;
      Exponent_Sign : Natural                    := 0;
      Base          : GMP.Base_Arg range 2 .. 36 := 10;
      Negative      : Boolean                    := False;
      Offset        : Natural                    := 0;
      First, Last   : Natural;
      Function_Name : constant String            := "To_MPFR_Float : ";
      Based_Literal : constant String            := "Invalid Based Literal";
   begin

      if Image (Image'First) = '-' then
         Negative := True;
         Offset   := 1;
      elsif Image (Image'First) = '+' then
         Offset := 1;
      else
         while Image (Image'First + Offset) = ' ' loop
            Offset := Offset + 1;
         end loop;
      end if;

      -- 1) find if it is a based literal
      First := Index (Image, "#");
      if First /= 0 then
         -- It is a based literal !
         Last := Index (Image (First + 1 .. Image'Last), "#");
         if Last = 0 then
            raise Data_Error with Function_Name &
                                  " couldn't find ending # in based literal";
         end if;
         Base :=
            GMP.Base_Arg'Value (Image (Image'First + Offset .. First - 1));
         if Last < Image'Last then
            if Image (Last + 1) = 'E' or else Image (Last + 1) = 'e' then
               Exponent_Sign := Last + 1;
            else
               raise Data_Error with Function_Name & Based_Literal;
            end if;
         end if;
         First := First + 1;
         Last  := Last - 1;

      else
         -- It is a decimal literal !
         First := Image'First + Offset;
         Last  := Image'Last;
         for i in First + 1 .. Image'Last loop
            if Image (i) = 'E' or else Image (i) = 'e' then
               Last          := i - 1;
               Exponent_Sign := Last + 1;
               exit;
            end if;
         end loop;
      end if;

      -- 2) suppress "_", check the digits are correct
      declare
         subtype Index_T is Integer range 0 .. Last - First + 1;
         C_Image  : String (1 .. Index_T'Last);
         C_Last   : Index_T;
         C_String : chars_ptr;
      begin
         C_Last := 0;
         for i in First .. Last loop
            if Image (i) /= '_' then
               C_Last           := C_Last + 1;
               C_Image (C_Last) := Image (i);
               if Character'Pos (C_Image (C_Last)) < Character'Pos ('0')
                 or else Character'Pos (C_Image (C_Last)) >=
                         Character'Pos ('0') + Base
               then
                  if Image (i) /= '.' and then Image (i) /= ',' then
                     raise Data_Error with Function_Name &
                                           "Found an invalid digit in literal";
                  end if;
               end if;
            end if;
         end loop;

         if Exponent_Sign /= 0 then
            C_String := New_String (C_Image (1 .. C_Last) &
                        "@" & Image (Exponent_Sign + 1 .. Image'Last));
         else
            C_String := New_String (C_Image (1 .. C_Last));
         end if;
         if Decimal_Digits = Default_Precision then
            mpfr_init_set_str (mpf, C_String, Base, default_rounding_mode);
         else
            mpfr_init2 (mpf, Prec_in_Bits (Decimal_Digits));
            mpfr_set_str (mpf, C_String, Base, default_rounding_mode);
         end if;
         Free (C_String);
      end;

      if Negative then
         mpfr_neg (mpf, mpf, default_rounding_mode);
      end if;
      return (Ada.Finalization.Controlled with mpf);
   exception
      when Data_Error => raise;
      when others     => raise Data_Error with Function_Name & "Invalid Literal";
   end To_MPFR_Float;

   -----------
   -- Image --
   -----------

   function Image
     (Item          : MPFR_Float;
      Base          : Base_Range := 10;
      Full_Mantissa : Boolean    := False)
      return          String
   is
      MPFR_Image : Strings.chars_ptr;
      expptr     : aliased GMP.mp_exp_t_ptr := new GMP.Specific.mp_exp_t'(0);
      -- aliased because it will be modified by the function call !!
      n_Digits : size_t;
   begin

      if Full_Mantissa then
         n_Digits := 0;
      else
         n_Digits := size_t (Prec_in_Digits (mpfr_get_prec (Item.Value)));
      end if;

      MPFR_Image :=
         mpfr_get_str
           (str      => Strings.Null_Ptr,
            expptr   => expptr,
            base     => GMP.Base_Arg (Base),
            n_digits => n_Digits,
            op       => Item.Value,
            rnd      => default_rounding_mode);
      declare
         Mantissa : String := Strings.Value (MPFR_Image);
         use GMP;

         function Exponent return String is
            Exp : String := Specific.mp_exp_t'Image (expptr.all - 1);
         begin
            if expptr.all - 1 /= 0 then
               if Exp (1) = ' ' then
                  return "E" & Exp (2 .. Exp'Last);  -- To remove the first space !!
               else
                  return "E" & Exp;
               end if;
            else
               return "";
            end if;
         end Exponent;

      begin
         mpfr_free_str (MPFR_Image);
         if Base = 10 then
            if mpfr_sgn (Item.Value) < 0 then
               return '-' & Mantissa (2) &
                      "." & Mantissa (3 .. Mantissa'Last) & Exponent;
            else
               return ' ' & Mantissa (1) &
                      "." & Mantissa (2 .. Mantissa'Last) & Exponent;
            end if;
         else
            if mpfr_sgn (Item.Value) < 0 then
               return '-' & Base_Range'Image (Base) &
                      '#' & Mantissa (2) &
                      "." & Mantissa (3 .. Mantissa'Last) &
                      '#' & Exponent;
            else
               return Base_Range'Image (Base) &
                      '#' & Mantissa (2) &
                      "." & Mantissa (3 .. Mantissa'Last) &
                      '#' & Exponent;
            end if;
         end if;
      end;
   end Image;

   ---------------------
   -- NaN, +Inf, -Inf --
   ---------------------

   function Not_A_Number return MPFR_Float is
      result : mpfr_t;
   begin
      mpfr_init (result); -- set Result to NaN.
      return (Ada.Finalization.Controlled with result);
   end Not_A_Number;

   function Plus_Infinity return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init (Result);
      mpfr_set_inf (Result, +1);
      return (Ada.Finalization.Controlled with Result);
   end Plus_Infinity;

   function Minus_Infinity return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init (Result);
      mpfr_set_inf (Result, -1);
      return (Ada.Finalization.Controlled with Result);
   end Minus_Infinity;

   -----------------
   -- Constant Pi --
   -----------------

   function Pi return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init (Result);
      mpfr_const_pi (Result, default_rounding_mode);
      return (Ada.Finalization.Controlled with Result);
   end Pi;

   -------------------------------
   -- "=", "<", "<=", ">", ">=" --
   -------------------------------

   function "=" (Left, Right : MPFR_Float) return Boolean is
      Result : int;
   begin
      Result := mpfr_equal_p (Left.Value, Right.Value);
      if mpfr_erangeflag_p > 0 then
         raise MPFR_Range_Error;
      end if;
      return (Result /= 0);
   end "=";

   function "<" (Left, Right : MPFR_Float) return Boolean is
      Result : int;
   begin
      Result := mpfr_less_p (Left.Value, Right.Value);
      if mpfr_erangeflag_p > 0 then
         raise MPFR_Range_Error;
      end if;
      return (Result /= 0);
   end "<";

   function "<=" (Left, Right : MPFR_Float) return Boolean is
      Result : int;
   begin
      Result := mpfr_lessequal_p (Left.Value, Right.Value);
      if mpfr_erangeflag_p > 0 then
         raise MPFR_Range_Error;
      end if;
      return (Result /= 0);
   end "<=";

   function ">" (Left, Right : MPFR_Float) return Boolean is
      Result : int;
   begin
      Result := mpfr_greater_p (Left.Value, Right.Value);
      if mpfr_erangeflag_p > 0 then
         raise MPFR_Range_Error;
      end if;
      return (Result /= 0);
   end ">";

   function ">=" (Left, Right : MPFR_Float) return Boolean is
      Result : int;
   begin
      Result := mpfr_greaterequal_p (Left.Value, Right.Value);
      if mpfr_erangeflag_p > 0 then
         raise MPFR_Range_Error;
      end if;
      return (Result /= 0);
   end ">=";

   ---------------------------------
   -- Initialize Result Precision --
   ---------------------------------

   function Initialize_Result_for_add
     (Left, Right : MPFR_Float)
      return        mpfr_t
   is
      Greatest, Smallest, Result     : mpfr_t;
      Precision, Smallests_Precision : mpfr_prec_t;
   begin
      if mpfr_cmpabs (Left.Value, Right.Value) > 0 then
         -- Greatest and Smallest are simple alias !
         Greatest := Left.Value;
         Smallest := Right.Value;
      else
         Greatest := Right.Value;
         Smallest := Left.Value;
      end if;
      Precision := mpfr_get_prec (Greatest);
      if mpfr_zero_p (Smallest) = 0 then
         Smallests_Precision := mpfr_get_prec (Smallest) +
                                mpfr_prec_t'Base (mpfr_get_exp (Greatest) -
                                                  mpfr_get_exp (Smallest));
         if Precision > Smallests_Precision then
            Precision := Smallests_Precision;
         end if;
      end if;
      mpfr_init2 (Result, Precision);
      return Result;
   end Initialize_Result_for_add;

   function Initialize_Result_for_mul
     (Left, Right : MPFR_Float)
      return        mpfr_t
   is
      Result                     : mpfr_t;
      Precision, Right_Precision : mpfr_prec_t;
   begin
      Precision       := mpfr_get_prec (Left.Value);
      Right_Precision := mpfr_get_prec (Right.Value);
      if Precision > Right_Precision then
         Precision := Right_Precision;
      end if;
      mpfr_init2 (Result, Precision);
      return Result;
   end Initialize_Result_for_mul;

   -------------------------
   -- Check for Exception --
   -------------------------

   procedure Check_For_Exception is
   -- This procedure must be used only inside MPFR.Floats body and not
   -- exported. It must be called just after the call of a MPFR function
   begin
      if mpfr_underflow_p   > 0 then raise MPFR_Underflow;
      elsif mpfr_overflow_p > 0 then raise MPFR_Overflow;
      elsif mpfr_nanflag_p  > 0 then raise MPFR_Not_a_Number;
      --elsif mpfr_inexflag_p > 0 then raise MPFR_Inexact;
      end if;
   end Check_For_Exception;

   ---------
   -- abs --
   ---------

   function "abs" (Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Right.Value));
      mpfr_abs (Result, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "abs";

   ---------
   -- "-" --
   ---------

   function "-" (Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Right.Value));
      mpfr_neg (Result, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "-";

   ---------
   -- "+" --
   ---------

   function "+" (Left, Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t := Initialize_Result_for_add (Left, Right);
   begin
      mpfr_add (Result, Left.Value, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "+";

   ---------
   -- "-" --
   ---------

   function "-" (Left, Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t := Initialize_Result_for_add (Left, Right);
   begin
      mpfr_sub (Result, Left.Value, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "-";

   ---------
   -- "*" --
   ---------

   function "*" (Left, Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t := Initialize_Result_for_mul (Left, Right);
   begin
      mpfr_mul (Result, Left.Value, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "*";

   ---------
   -- "/" --
   ---------

   function "/" (Left, Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t := Initialize_Result_for_mul (Left, Right);
   begin
      mpfr_div (Result, Left.Value, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "/";

   ---------
   -- "**" --
   ---------

   function "**"
     (Left  : MPFR_Float;
      Right : Integer'Base)
      return  MPFR_Float
   is
      Result    : mpfr_t;
      Precision : mpfr_prec_t := mpfr_get_prec (Left.Value);
      package Functions is new Ada.Numerics.Generic_Elementary_Functions (
         Float);
      use Functions;
   begin
      -- (x + e)**n = x**n + e*n*x**(n-1) + ...
      -- [ (x+e)**n - (x**n) ] / x**n = n * [(x+e) - x] / x
      -- => the relative error is multiplied by n, precision is reduced by log2(n) !
      if Right > 1 then
         if Right = 2 then
            Precision := Precision - 1;
         elsif Precision - mpfr_prec_t (Log (Float (abs (Right)), 2.0)) >
               mpfr_prec_t'First
         then
            Precision := Precision -
                         mpfr_prec_t (Log (Float (abs (Right)), 2.0));
         else
            Precision := mpfr_prec_t'First;
         end if;
      end if;
      mpfr_init2 (Result, Precision);
      mpfr_pow_si (Result, Left.Value, long (Right), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "**";

   ----------------------------
   -- With Standard Integers --
   ----------------------------

   function To_MPFR_Float
     (Value          : Integer'Base;
      Decimal_Digits : Natural := Default_Precision)
      return           MPFR_Float
   is
      result : mpfr_t;
   begin
      if Decimal_Digits = Default_Precision then
         mpfr_init (result);
      else
         mpfr_init2 (result, Prec_in_Bits (Decimal_Digits));
      end if;
      mpfr_set_si (result, long (Value), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with result);
   end To_MPFR_Float;

   ---------
   -- "+" --
   ---------

   function "+" (Left : MPFR_Float; Right : Integer'Base) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Left.Value));
      mpfr_add_si (Result, Left.Value, long (Right), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "+";

   function "+" (Left : Integer'Base; Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Right.Value));
      mpfr_add_si (Result, Right.Value, long (Left), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "+";

   ---------
   -- "-" --
   ---------
   function "-" (Left : MPFR_Float; Right : Integer'Base) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Left.Value));
      mpfr_sub_si (Result, Left.Value, long (Right), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "-";

   function "-" (Left : Integer'Base; Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Right.Value));
      mpfr_si_sub (Result, long (Left), Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "-";

   ---------
   -- "*" --
   ---------

   function "*" (Left : MPFR_Float; Right : Integer'Base) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Left.Value));
      mpfr_mul_si (Result, Left.Value, long (Right), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "*";

   function "*" (Left : Integer'Base; Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Right.Value));
      mpfr_mul_si (Result, Right.Value, long (Left), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "*";

   ---------
   -- "/" --
   ---------

   function "/" (Left : MPFR_Float; Right : Integer'Base) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Left.Value));
      mpfr_div_si (Result, Left.Value, long (Right), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "/";

   function "/" (Left : Integer'Base; Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Right.Value));
      mpfr_si_div (Result, long (Left), Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "/";

   --------------------------
   -- with Standard Floats --
   --------------------------

   function Generic_To_MPFR_Float (Value : F) return MPFR_Float is
   begin
      return To_MPFR_Float (F'Image (Value), F'Digits);
   end Generic_To_MPFR_Float;

   function To_MPFR_Float (Value          : Float;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float
   is
      result : mpfr_t;
   begin
      if Decimal_Digits = Default_Precision then
         mpfr_init (result);
      else
         mpfr_init2 (result, Prec_in_Bits (Decimal_Digits));
      end if;
      mpfr_set_flt (Result, Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with result);
   end To_MPFR_Float;


   function To_MPFR_Float (Value          : Long_Float;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float
   is
      result : mpfr_t;
   begin
      if Decimal_Digits = Default_Precision then
         mpfr_init (result);
      else
         mpfr_init2 (result, Prec_in_Bits (Decimal_Digits));
      end if;
      mpfr_set_d (Result, Double(Value), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with result);
   end To_MPFR_Float;

   -----------
   -- Round --
   -----------

   function Round (X : MPFR_Float) return Float
   is
      Result : Float;
   begin
      Result := mpfr_get_flt(X.Value, default_rounding_mode);
      return Result;
   end Round;

   function Round (X : MPFR_Float) return Long_Float
   is
      Result : Long_Float;
   begin
      Result := Long_Float(mpfr_get_d(X.Value, default_rounding_mode));
      return Result;
   end Round;

   function Round (X : MPFR_Float) return Long_Long_Float
   is
      Result : Long_Long_Float;
   begin
      Result := Long_Long_Float(mpfr_get_ld(X.Value, default_rounding_mode));
      return Result;
   end Round;

   function Generic_Round (X : MPFR_Float) return F
   is begin
      case F'Base'Digits is
      when Float'Digits           => return F(mpfr_get_flt(X.Value, default_rounding_mode));
      when Long_Float'Digits      => return F(mpfr_get_d  (X.Value, default_rounding_mode));
      when Long_Long_Float'Digits => return F(mpfr_get_ld (X.Value, default_rounding_mode));
      when others                 => raise Constraint_Error;
      end case;
   end Generic_Round;



   -------------
   -- Compare --
   -------------

   function Compare (Value : Float; Exact_Value : MPFR_Float) return Float is
      Difference : mpfr_t;
      Result     : Float;
   begin
      mpfr_init2 (Difference, mpfr_get_prec (Exact_Value.Value));
      mpfr_set_flt (Difference, Float (Value), default_rounding_mode);
      mpfr_sub
        (Difference,
         Difference,
         Exact_Value.Value,
         default_rounding_mode);
      Check_For_Exception;
      Result := Float (mpfr_get_flt (Difference, default_rounding_mode));
      mpfr_clear (Difference);
      return Result;
   end Compare;

   function Compare
     (Value       : Long_Float;
      Exact_Value : MPFR_Float)
      return        Long_Float
   is
      Difference : mpfr_t;
      Result     : Long_Float;
   -- We make the asumption that Long_Float == Double, which is true for GNAT !
   begin
      mpfr_init2 (Difference, mpfr_get_prec (Exact_Value.Value));
      mpfr_set_d (Difference, double (Value), default_rounding_mode);
      mpfr_sub
        (Difference,
         Difference,
         Exact_Value.Value,
         default_rounding_mode);
      Check_For_Exception;
      Result := Long_Float (mpfr_get_d (Difference, default_rounding_mode));
      mpfr_clear (Difference);
      return Result;
   end Compare;

   -----------------------------
   -- with Unbounded_Integers --
   -----------------------------

   function To_MPFR_Float
     (Value          : Unbounded_Integer;
      Decimal_Digits : Natural := Default_Precision)
      return           MPFR_Float
   is
      result : mpfr_t;
   begin
      if Decimal_Digits = Default_Precision then
         mpfr_init (result);
      elsif Decimal_Digits = Best_Precision then
         declare
            prec : mpfr_prec_t :=
               mpfr_prec_t (GMP.mpz_sizeinbase
                               (GMP.Integers.Aliased_Internal_Value (Value),
                                2) +
                            4);
         -- I need to add 4 bits to get the proper result, but i don't know why !
         begin
            if prec > mpfr_get_default_prec then
               mpfr_init2 (result, prec);
            else
               mpfr_init (result);
            end if;
         end;
      else
         mpfr_init2 (result, Prec_in_Bits (Decimal_Digits));
      end if;

      mpfr_set_z
        (result,
         GMP.Integers.Aliased_Internal_Value (Value),
         default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with result);
   end To_MPFR_Float;

   --------------------------------
   -- Round_To_Unbounded_Integer --
   --------------------------------

   function Round_To_Unbounded_Integer
     (Item : MPFR_Float)
      return Unbounded_Integer
   is
      mpz : GMP.mpz_t;
   begin
      GMP.mpz_init (mpz);
      mpfr_get_z (mpz, Item.Value, default_rounding_mode);
      return GMP.Integers.To_Unbounded_Integer (mpz);
   end Round_To_Unbounded_Integer;

   function Round_To_Unbounded_Integer
     (Item     : MPFR_Float;
      Rounding : Rounding_Mode_Kind)
      return     Unbounded_Integer
   is
      rounding_mode : mpfr_rnd_t := To_mpfr_rnd_t (Rounding);
      mpz           : GMP.mpz_t;
   begin
      GMP.mpz_init (mpz);
      mpfr_get_z (mpz, Item.Value, rounding_mode);
      return GMP.Integers.To_Unbounded_Integer (mpz);
   end Round_To_Unbounded_Integer;

   ----------------------------------------------------------------------------
   ---
   -- ELEMENTARY FUNCTIONS
   ----
   ----------------------------------------------------------------------------
   ---

   ----------
   -- Sqrt --
   ----------

   function Sqrt (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_sqrt (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Sqrt;

   ---------
   -- Log --
   ---------

   function Log (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_log (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Log;

   ---------
   -- Log --
   ---------

   function Log (X : MPFR_Float; Base : Integer'Base) return MPFR_Float is
      Result : mpfr_t;
      prec   : mpfr_prec_t := mpfr_get_prec (X.Value);
   begin
      mpfr_init2 (Result, prec);
      if Base = 2 then
         mpfr_log2 (Result, X.Value, default_rounding_mode);
      elsif Base = 10 then
         mpfr_log10 (Result, X.Value, default_rounding_mode);
      else
         declare
            Log_X, Log_Base, MPFR_Base : mpfr_t;
         begin
            mpfr_init2 (MPFR_Base, prec);
            mpfr_set_si (MPFR_Base, long (Base), default_rounding_mode);
            mpfr_init2 (Log_X, prec);
            mpfr_log (Log_X, X.Value, default_rounding_mode);
            mpfr_init2 (Log_Base, prec);
            mpfr_log (Log_Base, MPFR_Base, default_rounding_mode);
            mpfr_div (Result, Log_X, Log_Base, default_rounding_mode);
         end;
      end if;
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Log;

   ---------
   -- Exp --
   ---------

   function Exp (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_exp (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Exp;

   --------
   -- ** --
   --------

   function "**"
     (Left  : MPFR_Float;
      Right : Unbounded_Integer)
      return  MPFR_Float
   is
      Result   : mpfr_t;
      N        : Integer;
      Is_valid : Boolean;
   begin
      -- Since the exponant of a mpfr_t has a range similar to an integer,
      -- we simply convert the Unbounded_Integer into an integer...
      GMP.Integers.To_Standard_Integer (Right, N, Is_valid);
      if Is_valid then
         return Left ** N;
      else
         -- Do we need to keep this ??
         mpfr_init2 (Result, mpfr_get_prec (Left.Value));
         mpfr_pow_z
           (Result,
            Left.Value,
            GMP.Integers.Aliased_Internal_Value (Right),
            default_rounding_mode);
         Check_For_Exception;
         return (Ada.Finalization.Controlled with Result);
      end if;
   end "**";

   function "**" (Left, Right : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (Left.Value));
      mpfr_pow (Result, Left.Value, Right.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end "**";

   ---------
   -- Sin --
   ---------

   function Sin (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_sin (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Sin;

   ---------
   -- Cos --
   ---------

   function Cos (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_cos (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Cos;

   ---------
   -- Tan --
   ---------

   function Tan (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_tan (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Tan;

   ---------
   -- Sec --
   ---------

   function Sec (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_sec (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Sec;

   ---------
   -- Csc --
   ---------

   function Csc (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_csc (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Csc;

   ---------
   -- Cot --
   ---------

   function Cot (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_cot (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Cot;

   ------------
   -- Arcsin --
   ------------

   function ArcSin (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_asin (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end ArcSin;

   ------------
   -- Arccos --
   ------------

   function Arccos (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_acos (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Arccos;

   ------------
   -- Arctan --
   ------------

   function Arctan (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_atan (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Arctan;

   ----------
   -- Sinh --
   ----------

   function Sinh (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_sinh (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Sinh;

   ----------
   -- Cosh --
   ----------

   function Cosh (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_cosh (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Cosh;

   ----------
   -- Tanh --
   ----------

   function Tanh (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_tanh (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Tanh;

   ----------
   -- Sech --
   ----------

   function Sech (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_sech (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Sech;

   ----------
   -- Csch --
   ----------

   function Csch (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_csch (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Csch;

   ----------
   -- Coth --
   ----------

   function Coth (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_coth (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Coth;

   -------------
   -- Arcsinh --
   -------------

   function ArcSinh (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_asinh (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end ArcSinh;

   -------------
   -- Arccosh --
   -------------

   function Arccosh (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_acosh (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Arccosh;

   -------------
   -- Arctanh --
   -------------

   function Arctanh (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_atanh (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Arctanh;

   ---------------
   -- Factorial --
   ---------------

   function Factorial (N : Natural'Base) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init (Result);
      mpfr_fac_ui (Result, unsigned_long (N), default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Factorial;

   -----------
   -- Gamma --
   -----------

   function Gamma (X : MPFR_Float) return MPFR_Float is
      Result : mpfr_t;
   begin
      mpfr_init2 (Result, mpfr_get_prec (X.Value));
      mpfr_gamma (Result, X.Value, default_rounding_mode);
      Check_For_Exception;
      return (Ada.Finalization.Controlled with Result);
   end Gamma;

   ----------------
   -- Initialize --
   ----------------

   procedure Initialize (Object : in out MPFR_Float) is
   begin
      mpfr_init (Object.Value);
   end Initialize;

   ------------
   -- Adjust --
   ------------

   procedure Adjust (Object : in out MPFR_Float) is
      mpf : mpfr_t := Object.Value;
   begin
      mpfr_init2 (Object.Value, mpfr_get_prec (mpf));
      mpfr_set (Object.Value, mpf, default_rounding_mode);
      if required_min_precision > mpfr_prec_t'First
        and then required_min_precision < mpfr_get_prec (mpf)
      then
         raise Not_Accurate;
      end if;
   end Adjust;

   --------------
   -- Finalize --
   --------------

   procedure Finalize (Object : in out MPFR_Float) is
   begin
      mpfr_clear (Object.Value);
   end Finalize;

end MPFR.Floats;
